Clustering of gyrotactic microorganisms in turbulent flows 
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We study the spatial distribution of gyrotactic microorganisms transported by a three-dimensional 
turbulent flow generated by direct numerical simulations. We find that gyrotaxis combines with 
turbulent fluctuations to produce small scales (multi-)fractal clustering. We explain this result 
by showing that gyrotactic swimming cells behave like tracers in a compressible flow. The effective 
compressibility is derived in the limits of fluid acceleration much larger and smaller than the gravity. 

PACS numbers: 05.45.-a, 47.63. Gd, 92.20.jf 



Microbial patchiness in oceans is important for ecologi- 
cal and evolutionary dynamics [HIH and for biogeochemi- 
cal processes Q. In motile aquatic microorganisms, self- 
propulsion provides a mechanism to escape fluid path- 
lines, potentially leading to small-scale patchiness [iHSl- 
Remarkably, motility combined with fluid flows can also 
generate large-scale inhomogeneities. For instance, spec- 
tacular aggregation of phytoplankton cells (in layers cen- 
timeters to meters thin, horizontally extending from hun- 
dreds of meters to kilometers) can result from vertical 
shears and gyrotactic swimming Q- Gyrotaxis charac- 
terizes several species of motile microalgae whose swim- 
ming direction is determined by the balance of viscous 
and gravitational torques, due to the displacement be- 
tween the cell center of mass and buoyancy. As an effect 
of such balance, for example, gyrotactic algae aggregate 
in the center (wall) of descending (ascending) vertical 
pipe flows [1, H|. Gyrotaxis is observed in algae, e.g., of 
the genus Chlamydomonas, which can be engineered to 
trans por t microloads , or Dunaliella, employed in bio- 
fuels So far most studies focused on the dynamics 
of gyrotactic microorganisms in simple stationary flows 




or kinematic models [6H9|, 12, 13j | . 

In this Letter, we investigate the interplay between 
gyrotactic motility and realistic turbulent flows, as oc- 
curring in the sea. We find that turbulence and gyro- 
taxis combine to generate inhomogeneous distributions 
with small-scale (multi-)fractal statistics (see Fig. [TJ. 
We study the limit of gravitational acceleration much 
smaller or larger than turbulent accelerations to identify 
the mechanisms responsible for gyrotactic clustering in 
terms of an effective compressible velocity field. 

We consider dilute suspensions of non interacting 
motile microorganisms, much smaller than the smallest 
scale of turbulence, the Kolmogorov length n. We can 
thus model them as self-propelled particles with velocity, 



FIG. 1. (color online) Spatial distribution of gyrotactic 
swimmers (dots) in a slab of a 3D turbulent flow. Color 
code: yellow/blue corresponds to high/low vorticity values 
(In |w|/o; rm s). (Left) Limit of orientation dominated by local 
fluid acceleration (A — a, see text) with aggregation in high 
vorticity regions. (Right) Limit of gravity dominated orien- 
tation (A = —g). Parameters correspond to circled symbols 
in Fig. [2ji and Fig. [3^l, respectively. 



are assumed spherical and neutrally buoyant, with the 
center of mass displaced by h with respect to the geo- 
metric one. The swimming direction p, determined by 
the total torque acting on the cell, evolves as 
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[A- (A-p)p] + -u> x p, 



(2) 



where u> is the fluid vorticity and v — Zv /h is the orien- 
tation speed for spherical cells subject to the acceleration 
j4.[9|]. In a fluid at rest, besides viscous forces, only grav- 
ity (and buoyancy) g is acting and thus A = —g = gz, 
while acceleration due to swimming is neglected In 
presence of a flow, we have A = a — g where 



dtu + u ■ Vu = — Vp + vV 2 u + / 



(3) 



X = U(X,t)+V S P; 



(1) 



given by the sum of the fluid velocity u at the particle 
position X and the swimming contribution VsV, where 
the swimming speed v s is assumed constant [a, |9| • Cells 



is the fluid acceleration given by the Navier-Stokes equa- 
tions ruling the velocity u of an incompressible (V • u = 
0) fluid with viscosity v, pressure p and stirred by an 
external forcing /. Previous studies on gyrotactic swim- 
mers disregarded fluid acceleration, as mainly focused 
on simple, non-turbulent flows where \a\ <C g. In tur- 
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bulcnce, fluid acceleration can locally exceed g 14| and 
therefore its contribution has to be taken into account. 

The first term on the rhs of Eq. @ causes the di- 
rection of swimming p to align with A on a time scale 
v l A. When the contribution of fluid acceleration can 
be neglected, cells tend to orient vertically (p — > z) on 
a time scale B = v a /g. The alignment is contrasted by 
the vorticity term uj x p and, depending on Bui being 
smaller or larger than 1, cells may swim along a result- 
ing local equilibrium direction or tumble randoml y as the 
orientation becomes unstable due to vorticity P. 1l3j]. In 
principle, the swimming direction may be modified also 
by rotational Brownian motion [l5| and tumb ling due to 
flagella desynchronization during swimming [16| , which 
are here neglected. The former effect is very small for 
typical algae (having size O(10/xm)); the latter can be 
neglected whenever the tumbling time is longer than the 
reorientation one. 

We study gyrotactic swimming in homogeneous and 
isotropic turbulent velocity fields of moderate intensity 
(Rex ~ 65 — 100) by means of direct numerical simula- 
tions of Navier-Stokes equations. In particular, Eq. Q 
is solved by means of a standard pscudospcctral algo- 
rithm with 2 nd order Runge-Kutta time-stepping, on a 
tri-periodic cubic grid of size N 3 (for AT = 128 and 256). 
Statistical stationarity is guaranteed by means of a zero 
mean, Gaussian and white in time random forcing / re- 
stricted to large scales. Viscosity v is such that the Kol- 
mogorov length r\ is of the order of the grid spacing, en- 
suring well resolved small-scale velocity dynamics. For 
different values of g, several populations of swimmers, 
characterized by different values of v s and v are injected 
with random positions and orientations. At each time 
step, velocity and acceleration at the swimmers positions, 
needed to integrate Eqs. ([T][2]), are obtained by inter- 
polation. The self-propelled particles are then evolved, 
and their distribution and orientation studied in statis- 
tically steady conditions. In the sequel, we mostly focus 
on the dependence on the orientation speed v by fixing 
v s rj 0.3u r; , u v being the typical fluid velocity fluctuation 
at the Kolmogorov scale. 

Formally, Eqs. (Q~|[2]) define a dissipativc dynamical sys- 
tem evolving in the 2d-dimensional (actually 2d — 1 be- 
cause p 2 = 1 and d = 3) phase space (X,p) with phase- 
space contraction rate 
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As p orients in the direction a — g, T is expected to be 
negative on average, meaning that swimmers will evolve 
onto a dynamical attractor of dimension smaller than the 
whole phase space, which explains why clustering can 
be observed: if the fractal dimension of the attractor is 
smaller than d, clustering in position space (as in Fig. Q]) 
is possible (see Ref. [I?} for a conceptually similar phe- 
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FIG. 2. Swimmer properties as a function of the orientation 
parameter Bu> Ims (B — v /a Tms ) in the limit \a\ 3> g- (a) 
Average alignment with fluid acceleration (6 • p). (b) Aver- 
age square vorticity at swimmers position normalized to the 
volume average value, (c) correlation dimension Di. Circled 
symbol in (c) corresponds to the distribution shown in Fig.[T^i. 



nomenon occurring for inertial particles). We remark 
that clustering is a consequence of swimming: indeed 
for v s — Eqs. JT]) and @ decouple, thus cells become 
tracers advected by an incompressible velocity and can- 
not cluster. Moreover, in the limit v Q — y oo we have 
r — y and therefore swimmers cannot cluster. Nonethe- 
less, even in this limit, if v s > they deviate from fluid 
trajectories and generate interesting dynamics [l8j |. 

We now discuss the physical mechanisms of clustering 
which, as anticipated, depend on whether the dominating 
effect comes from the gravitational (g) or fluid accelera- 
tion (which we quantify in terms of its rms value a rms ). 

We start considering the case a rms ^> g and therefore 
we take A = a in Eq. J2]). Figure [2] summarizes the be- 
havior of the main obscrvables as a function of the dimen- 
sionless number j5w rms (now B = v /a rms ) measuring the 
ratio of the alignment timcscale to rotation timescale in- 
duced by vorticity. When the alignment is very fast, the 
swimming direction p becomes parallel to the local di- 
rection of the fluid acceleration a = a/ a, as confirmed 
by Fig. [5h. showing that (d ■ p) — y 1 for Bui Tms <C 1 (here 
and in the following ([•]) denotes average over particle 
distribution). In this limit, swimming cells behave like 
tracers advected by an effective velocity v m u + v s a. 
While u is incompressible, the effective velocity field v 
is not: V • v oc v s V • a being negative (positive) in high 
vorticity (strain) regions. Therefore, as it occurs for in- 
ertial particles lighter than fluid [l9|, H3] , the swimmers 
cluster inside vortical structures (Fig. [Tk. and Fig. [5b). 
The divergence of v is proportional to v s , clustering is 
thus expected to increase with the swimming speed. In 
the opposite limit of slow alignment, when Bu! rms ^> 1, 
random tumbling due to fluid vorticity dominates, hence 
swimming orientation cannot align to the local accelera- 
tion ((a ■ p) —y0, sec Fig. [5^): the compressible effect is 
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lost and particles distribute uniformly in the volume. To 
quantify clustering we measured the correlation dimen- 
sion D 2l ruling the small-distance (r — > 0) behavior of 
the probability to find two swimmers at separation less 
than r: P 2 (\X 1 - X 2 \ < r) oc r° 2 [2l|. For uniformly 
distributed particles D 2 = d, while when clustering is 
present the probability to find close pairs increases and 
D 2 < d (see e.g. [22| for a similar study in the case of 
inertial particles). In Fig.[2j: we show D 2 as a function of 
Bu> rms : for Bui lms <C 1, D 2 w 1.5, indicating strong clus- 
tering in almost filamental structures; conversely, when 
Buj rms > 1, the correlation dimension approaches the 
uniform-distribution value D 2 w 3. 

We now consider the limit a lms -C g when we can take 
A = g and Eq. © reads 

1 s 1 

P = 7^( z -P*PJ + x P« ( 5 ) 

with B = v j g. Similarly to the previous case, when 
Buj rms — ► the cells orient in the preferred direction z, 
which is now fixed in space. The effective velocity thus 
becomes v = u + v s z which, unlike the previous case, is 
incompressible (V • v = 0). Therefore, now we expect 
that not only for Buj Tms 3> 1 but also for Bu> lms — > 
swimmers distribute uniformly, as confirmed by Fig. [3^ 
showing that D 2 — > 3 in both limits. Remarkably, Fig.[3K 
shows that also in this case gyrotactic swimmers clus- 
ter on a fractal set (see Fig. [Tb) for intermediate values, 
with a well defined minimum of the correlation dimension 
(D 2 « 2.7) for Buj lms ~ 0(1). We remark than an opti- 
mal orientation timescale for aggregation is also observed 
in steady kinematic vortical flows [6| where, however, a 
vast class of trajectories is integrablc. 

We can understand the origin of the observed cluster- 
ing by considering the limit Buj lms < 1. In such limit, 
cell orientation being very fast we can assume that the 
swimming direction p is always at an equilibrium orien- 
tation with p x ,p y <C Pz — 1 (see Fig. [3Jd). In particular, 




FIG. 3. (color online) Clustering properties as a function of 
Bwrms, for g ^> \a\ (B — vo/g). (a) Correlation dimension 
D2 of the swimmer positions. Circled symbol corresponds to 
the data shown in Fig. [Tb. (b) Variances of swimming direc- 
tion components ((pj) = (p y ) in red, and (p^)). The dashed 
blue curve is the parabola (Buj TlnB ) 2 . The solid horizontal line 
represents the random orientation value 1/3. 
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FIG. 4. (color online) (a) Average fluid acceleration (positive, 
red filled circles) and velocity (negative, blue open circles) 
along the vertical at swimmer positions, expressed in per- 
centage of a r ms and v s , respectively, (b) Correlation — {p y u!x) 
vs Bw rms . The maximum for Ba; rms ~ 0(1) is understood 
noticing that — (p v cj m ) must decrease for -Bcj rms 3> 1. In 
the limit Bw rms C 1, —(p y uj x ) ~ -B^r ms (solid line) is im- 
plied by p ~ (Boj x , —Bujy, 1), see text. In the random 
tumbling limit (p^) = (p y ) — > 1/3 (Fig. [3^), which implies 
— {pyUJ x ) = (j>y) 3 /B ~ 1/(3B) (dashed line). 

solving Eq. ([5]) for p = 0, at first order in p Xl p y , one 
finds p x ~ Buiy and p y ~ —Buj x (which is confirmed by 
simulations). As a consequence, the effective swimmer 
velocity field v = u + v s p with p ~ {Bu> y , —Bu> x , 1) has 
a compressible component with divergence 

V • v ~ ~v s BV 2 u z , (6) 

which, unlike the previous case, is unrelated to fluid ac- 
celeration so that swimmers will cluster in regions dif- 
ferent from those of high vorticity (compare Fig. [T^, and 
b). We notice that ([6]) generalizes the well known mech- 
anism of cell focusing in the center (walls) of downward 
(upward) vertical pipe flows Q . Notice that in the above 
argument the vertical component of the vorticity plays 
no role, as it does not change p z . 

Another consequence of the expansion p ~ 
(Bio x ,—B(jj y ,l) is that p x (resp. p y ) and lo v (lo x ) 
have locally the same (opposite) sign. Numerical 
simulations show that this remains true also for larger 
values of Buj rms , on average. Indeed, at stationarity, by 
averaging Eq. ([5]) and using isotropy on the (x, y) plane 
(guaranteed by the isotropy of the fluid velocity field) 
we obtain (p£) = {p 2 y ) = B{p x uj y ) = -B(p y u> x ). The 
correlation between the horizontal components of p and 
us implies that the swimmers will stay longer in regions 
of the flow characterized by positive vertical velocity 
and negative vertical acceleration (Fig.QJi). This can be 
easily seen in a case with, say, a vortex aligned with the 
x-axis, where the above argument with u> x > implies 
(pz) > 0, (p y ) < 0, so that the trajectories spend more 
time in regions where a z > 0,u z < as there the swim- 
ming velocity opposes that of the fluid. The preferential 
concentration in these regions of the flow will be maximal 
(and correspondingly the correlation dimension minimal, 
i.e. clustering stronger) for Bu lms ~ 0(1) where the 
correlation between swimming direction and vorticity 
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FIG. 5. (color online) (a) Average square vorticity at swim- 
mers position normalized to the volume average value at vary- 
ing the ratio a = drms/s (a = corresponds to data of Fig.[3|; 
(b) correlation dimension D2 vs a; (c) D2 vs swimming speed 
v 3 /u v , the circled symbol corresponds to the circled one in 
(b); (d) generalized dimensions D q vs q for circled data in 
(b), notice that the case a — (filled black circles) appears 
to be less multifractal than when also the fluid acceleration is 
contributing to clustering (empty blue circles) . 



— {pyLo x ) = (p x io y ) is also maximal (Fig. Hb). For such 
value of Buj rms fluid regions with maximal deviation of 
the swimming direction from the vertical will balance 
vorticity dominated ones where p tumbles randomly. 
We observe that the swimmer vertical migration can 
be strongly inhibited by the bias towards downwelling 
regions: in Fig. [4ji, e.g., (u z ) can reach 30% of the 
swimming speed v s . 

In the general case, the relative importance of fluid 
and gravitational accelerations for clustering depends on 
the ratio a = a rms / g. Figure [5b indeed shows that the 
bias towards regions of high vorticity decreases with a 
and is absent when only the gravitational torque is act- 
ing (a = 0). The correlation dimension D2, shown in 
Fig. [5p, smoothly varies with a, interpolating from the 
two limits shown in Fig. [2k and [3b. We observe that, as 
anticipated, clustering is more effective for large swim- 
ming speeds as displayed in Fig. [5b, showing that, at 
fixed value of Bu rias , D2 decreases with v s /u n . Finally, 
as one can expect from general considerations on dynam- 
ical attractors [l?} , Fig. [5ji demonstrates that the spatial 
distribution of the gyrotactic self-propelled particles is 
multifractal, as the generalized dimensions D q (control- 
ling the probability to find q pa rticles at small separation) 



depends on the moment q [21 



identified two mechanisms driving microorganism clus- 
tering: the focusing in vortical regions due to local ad- 
justment of the swimming orientation with fluid accelera- 
tion, and the correlation between vorticity and swimming 
direction on the plane perpendicular to gravity leading 
particles to preferentially explore downwelling, upward 
accelerating regions. In general, gravity is expected to 
dominate when turbulent intensity is not very high and 
it is likely the most important effect in the ocean. Cru- 
cial parameters for observing clustering are in this case 
the ratio between swimming speed and small-scale fluid 
velocity fluctuations (v s /u v ) and the reorientation time 
scale with respect to vorticity intensity (-Bw rms ). For 
typical microalgae B « 1 — 6s and v s = 100 — 200^m/s 

[Ennui- In the 

ocean, the turbulence intensity, mea- 



sured in terms of kinetic energy dissipation e, varies 
from e ~ 10~ — lO~ 5 W/Kg in the upper mixing layer 
down to e ~ 10~ 6 — 10~ 7 W/Kg a few meters deeper 
M [ii| . We can thus estimate that V s /u n £ [0.02 : 0.4] 
and Buj rms £ [0.1 : 50] therefore the effects discussed in 
this Letter are relevant in realistic conditions and can 
definitely be tested in laboratory by tuning turbulence 
characteristics. 

We conclude by remarking that for non-spherical cells 
such as, e.g., prolate spheroids the term 7p-§ - (I — p®p) 
should be added to Eq. ([2]) (7 being the eccentricity, and 
S and I the symmetric rate of strain tensor and identity 
matrix, resp.) Such term is also contributing to the 
phase-space contraction rate (j?]) providing an additional 
mechanism for clustering [f|. It will thus be interesting 
to study if and how gyrotactic clustering in turbulence is 
modified at varying the cell shape. 
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Summarizing, wc have shown that gyrotactic motil- 
ity and realistic turbulent flows can generate small-scale 
patchiness (down to the Kolmogorov scale) in the distri- 
bution of bottom-heavy swimming microorganisms. Wc 
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